Scale invariance in early embryonic development

The body plan of the fruit fly is determined by the expression of just a handful of genes. We show that the spatial patterns of expression for several of these genes scale precisely with the size of the embryo. Concretely, discrete positional markers such as the peaks in striped patterns have absolute positions along the anterior-posterior axis that are proportional to embryo length, with better than 1% accuracy. Further, the information (in bits) that graded patterns of expression provide about position can be decomposed into information about fractional or scaled position and information about absolute position or embryo length; all of the available information is about scaled position, again with ~ 1% accuracy. These observations suggest that the underlying genetic network exhibits scale invariance in a deeper mathematical sense. Taking this mathematical statement seriously requires that the network dynamics have a zero mode, which connects to many other observations on this system.


I. INTRODUCTION
Closely related organisms can vary widely in size, but variations in their proportions are much smaller [1][2][3].There is a considerable gap between this qualitative observation and some precise mathematical statement of scaling, e.g., that the linear dimensions of all elements in the body plan are in direct proportion to the linear dimensions of the organism.If correct this scale invariance would be visible not only in the fully developed organism but already at some earlier stages in its development.
There are many examples of "allometric scaling," power-law relationships among different quantities across a well-defined class of organisms [4][5][6].In some cases, these relations connect the linear dimensions of different body parts.Nonetheless, truly precise spatial scaling in embryonic development would be quite surprising.
We understand the mechanisms of pattern formation in a wide range of non-biological systems, from fluid flows to crystal growth (snowflakes) and more [7][8][9][10][11][12], but none of these examples exhibit scale invariance.Instead, the elements of the pattern have linear dimensions set by microscopic parameters, and larger systems exhibit more repetitions of the same pattern rather than expansion or contraction of pattern elements to match the size of the system as a whole [13].Going back to the pioneering work of Turing [14], the mathematical structure of the equations governing these systems is not so different from the structure of models for genetic or biochemical networks.If we take these analogies literally, we would predict that taller people should have more vertebrae, which is obviously wrong.Is there a real problem here, or are we just oversimplifying?
Here we try to make the notion of scale invariance in development more precise.We use the first hours of development in the fruit fly as an example, following spatial patterns of morphogen concentration as they flow through three layers of a genetic network, from maternal inputs to the gap genes to the pair rule genes [15][16][17].In the spirit of earlier work [18][19][20][21][22] we analyze discrete positional markers, such as the stripes in pair rule-gene expression, and find that positions of these markers vary in proportion to the length of the embryo with better than 1% accuracy [23].We then go beyond discrete markers, decomposing the information carried by graded patterns of gap gene expression into information about fractional or scaled position vs. information about the absolute position; we find that all the available information is about fractional position along the anterior-posterior axis.Information that would signal a deviation from scale invariance is less than 1% of the total.
These results provide strong evidence for scaling in a precise mathematical sense, for both the gap genes and the pair rule genes.But at least one of the maternal inputs, Bicoid [24,25], does not show any sign of scale invariance: as in well-understood non-biological patternforming systems, there is a length scale that presumably is set by underlying molecular parameters and does not adjust in response to the linear dimensions of the whole embryo.This suggests that scale invariance is an emergent property of the gap gene network.
We argue that true scale invariance places very specific requirements on the dynamics of this network, independent of molecular details: it must have a "zero mode."This has connections to other observations on gap gene dynamics [26,27] and to more detailed models [28,29].

II. TESTING FOR SCALING
To make the notion of scaling more precise we take seriously the idea that cell fates are determined by the concentration of particular molecules called morphogens [30].Since the cell fates are tied to their positions, the concentrations of morphogens must also carry information about position along the body axes.These ideas are especially crisp in the early fly embryo, where we know the identities of all the relevant morphogens and rich spatial patterns in the concentrations of these molecules are established before cells make large-scale movements [31].
We focus on pattern formation along a single axis, which will be the anterior-posterior axis in the analysis of fly embryos below.Then we can measure position along this axis by a single variable 0 < x < L, where x = 0 is the anterior end of the embryo, x = L is the posterior end, and hence L is the length of the embryo.There are multiple morphogen species, indexed by i, and if we neglect the discreteness of cells then their concentration profiles are described by continuous functions g i (x; L).The notation emphasizes that concentration profiles may be different in embryos of different size L.
True scale invariance is the statement that the concentration of morphogens depends only on position relative to the length, that is If there is a fixed map from morphogen concentrations to cell fates, this scaling behavior would guarantee that cells adopt a fate that depends on their relative position x/L, and not separately on x and L. How do we test for scale invariance?If the concentration of the morphogen has a single peak as a function of x, we can write then scale invariance as in Eq. ( 1) requires that all the L dependence is contained in the position of the peak.
where f p is the fractional or scaled peak position, ⟨• • • ⟩ is the average over many embryos with different lengths, and noise allows that positions jitter from embryo to embryo.We emphasize that Eq. ( 3) is not just the statement that positional markers adjust (in absolute distance) to the length of the embryo; scale invariance as we have defined it in Eq. (1) requires that this adjustment is exactly linear with zero intercept.There is a natural generalization to concentration profiles that have multiple peaks, as with the pair rule genes (Fig. 1A, B).It has been known for some time that the morphogens in the early fly embryo carry enough information to specify scaled positions with ∼ 1% precision all along the anterior-posterior axis [32,33].At the same time, embryos from the same mother, in an inbred laboratory stock, fluctuate in length with a standard deviation of σ L /⟨L⟩ ∼ 4% (Appendix A and [34,35]).It would seem that to make these numbers consistent with one another, positional signals must scale with embryo length, but this is a bit subtle.
Imagine a hypothetical embryo in which, e.g., the peak of the morphogen profile is perfectly anchored in absolute position relative to the anterior pole of the embryo, with no scaling and no noise, such that x p = ⟨x p ⟩. Then the relative or fractional positions f p = x p /L fluctuate only because the lengths of the embryos vary, Thus for a marker that on average is a quarter of the way from the anterior to posterior, ⟨x p ⟩ = 0.25⟨L⟩, fluctuations will be σ fp (A) ∼ 0.01 even without scaling.Similarly, if we have a marker anchored at some fixed absolute position relative to the posterior then the variance in fractional position will be We can imagine cells combining anterior and posterior signals to reduce the error, With σ L /⟨L⟩ ∼ 0.04, fluctuations in fractional position thus could be less than ∼ 1.4% everywhere along the anterior-posterior axis, even in the absence of any scaling mechanism.Convincing ourselves that pattern formation is truly scale invariant requires a very precise measurement and depends on the system itself being very precise.
It is intuitive to think about scaling as the proportionality of positions to embryo length, as in Eq. ( 3), but it should be possible to test the scaling of the entire morphogen profile, as in Eq. ( 1), more directly.There are two related observations.First, to compare morphogen profiles across embryos of different lengths, we need a metric.Second, since morphogen profiles are noisy, it is unrealistic to expect the exact equality of two functions across all values of x.Fortunately, the noise level itself provides a metric for comparison, which is made precise in the language of information theory.
The statement that morphogen profiles depend on x and L means that the concentrations of these molecules provide information about the underlying positional variables.This is quantified, uniquely, by the Shannon information [36][37][38], I(g → {x, L}) = dg dx dL P (g|{x; L}) P (x, L) log 2 P (g|{x, L}) P (g) bits, where for more compact notation we write g = {g i } and dg = i dg i .Here P (g|{x, L}) is the probability of finding the set of morphogen concentrations {g i } at position x in an embryo of length L; P (g) is the probability of finding these concentrations averaged over all values of x and L; and P (x, L) is the distribution of positions and lengths.It is useful to recall that this information is mutual: the concentrations of morphogens provide cells with information about position, and specifying position allows us to predict the concentrations, so we write I(g; {x, L}).Information depends on both the mean spatial profiles of the morphogens and their noise levels.
True scale invariance would mean that all of the information conveyed by the morphogens is about the fractional position x/L: Equivalently, if we want to predict the morphogen concentration, it is enough to specify the fractional position, and no extra information is gained by knowing x and L separately.We can think of the total information as having a component about the relative position and an extra increment that describes the deviation from scaling, and we will see that with samples from a sufficiently large number of embryos, we can make a reliable estimate of ∆I.The smaller the fraction ∆I/I(g; x/L) the closer the system is to a mathematical ideal of scaling.More explicit expressions for ∆I are developed in Appendix B and applied to experiments in §IV.
We emphasize that true scale invariance, corresponding to ∆I = 0, is a very strong condition.Different levels of evidence for scaling in embryonic development have inspired models in which competing mechanisms can provide some cancellation of the intrinsic length scales determined by diffusion constants and reaction rates [39][40][41][42].These models typically allow for scaling in the position of a single discrete positional marker (e.g., the middle of the embryo), or for approximate scaling across a larger segment of the relevant axes.True scale invariance would require new dynamical mechanisms.

III. STRIPES AND BOUNDARIES
In the early fly embryo, information about position along the anterior-posterior axis flows from maternal inputs through the network of gap genes to the pair-rule genes [17].The pair-rule genes are expressed in interdigitating striped patterns that provide a preview of the segmented body plan in the fully developed organism; these stripes are visible within three hours after the egg is laid (Fig. 1A-C).The positions of pair-rule stripes are a clear example of the positional markers discussed above.
Here we analyze the spatial profiles of gene expression for three of the pair-rule genes-eve, prd, and runmeasured using fluorescent antibody staining of the corresponding proteins in more than one hundred embryos that were fixed during nuclear cycle 14 (nc14), i.e. between 2 and 3 h of development [33].Our results recapitulate earlier work [23] on a larger ensemble of embryos.
As soon as the stripes are visible it is straightforward to measure their positions x i [33].The time during nc14 can be measured with ∼ 1 min precision by following the progression of the invaginating cellularization membrane [32].The stripe positions vary systematically in time [19,[43][44][45][46] and are well described by as shown for the Eve stripes in Fig 1D .Combining data from all time points, we shift each embryo to the reference time t 0 = 45 min, We use this same procedure for the Prd and Run stripes, although these become clear only at slightly later times.
Figure 1E shows that the stripe positions x i measured from the anterior pole are proportional to the length of the embryo L.More precisely, if we fit these linear relations then intercepts are zero and slopes are equal to the mean fractional positions, as in Eq. (3), both results with error bars of less than 1% (Appendix C).This provides prima facie evidence for scaling of the pair-rule stripes, reinforcing the conclusions of earlier work [18][19][20][21].
We can go beyond the mean behaviors to look at fluctuations around these means.For each stripe i in each embryo α, we can write where ⟨• • • ⟩ now is an average over all the embryos in our sample.The variance of the relative position is σ 2 fi = ⟨(δf i ) 2 ⟩, and Fig. 1F shows that σ fi ≤ 0.01 for all 21 pair rule stripes that we measure.This is consistent with previous measurements, and with the information content of the gap gene expression patterns that feed into the generation of pair-rule stripes [33,47], but earlier work did not address scaling explicitly.
As a caution, we note that the observation of scaling in fixed embryos would be trivial if variations in embryo length were dominated by shrinkage during fixation.Across N em = 609 fixed embryos used for the analysis of gap genes (below) we find a mean length ⟨L⟩ fix = 455 µm, while across N em = 610 live embryos ( §V) we find ⟨L⟩ live = 490 µm.Hence, shrinkage with fixation is a bit less than 10% across many different experiments.But the variations in length are almost the same, (σ L /⟨L⟩) fix = 0.038, while (σ L /⟨L⟩) live = 0.037.The small extra variance in the length of fixed embryos cannot explain the scaling behavior that we observe.Figure 1F also shows that the fluctuations in fractional position are smaller than the bound on mechanisms that have no explicit scaling, from Eq. ( 7).This bound is very tight, because of the small variance in emrbyo lengths, and thus requires extreme precision in the measurement and biological reproducibility of the fractional positions to demonstrate scaling.To emphasize the importance of precision, we note that the position of the cephalic furrow is directly regulated by pair rule gene expression [48], but it has a slightly higher relative positional variance, due to the experimental difficulty of defining morphological features to less than the width of a single cell [49].
Here we show explicitly that the furrow position scales with embryo length (Appendix C).Even though the precision of the CF position is almost ∼ 1% in the scaled coordinates [49], this alone is not sufficient to reject the hypothesis that positions are defined in absolute rather than relative coordinates, as can be seen from Fig. 1F.
The pair rule stripes are shaped by input from the gap genes [50], and it is natural to ask whether the scaling behavior that we observe is inherited from these inputs.The gap genes were long discussed in terms of "expression domains," as if they were on/off switches [51][52][53][54].We now know that this misses a substantial fraction of the positional information encoded by these genes [33,47,55], but defining the boundaries of the expression domains as positional markers (Fig. 2A-D) allows us to give a preliminary analysis of scaling by following the same ideas as for the positions of the pair-rule stripes.
Previous experiments have measured the expression profiles of the gap genes [33], staining N em = 609 fixed embryos in nc14 with fluorescent antibodies directed at the proteins encoded by the gap genes (Fig. 2A-D).We define expression boundaries as the positions where the concentrations are half their maximum mean value, and we correct their relative positions to t 0 = 45 min as above.Figure 2E shows that all thirteen of the gap gene boundaries defined in this way have absolute positions that scale precisely with embryo length, as with the positions of the pair rule stripes.The accuracy of this scaling again is better than ∼ 1%, and this precision is better than the limiting performance of mechanisms that do not have some explicit sensitivity to embryo length (Fig. 2F).For the gap genes, this procedure allows us to span almost the entire range of the anterior-posterior axis.
In summary, stripes and boundaries of gene expression in the early fly embryo provide discrete positional markers, and the absolute positions of these markers are in all cases proportional to the length of the embryo.This is consistent with previous observations [18][19][20][21], but the precision of the scaling that we observe here is surprising.This suggests that the underlying genetic network exhibits true scale invariance, which we now test using the information decomposition [Eq.(10)].

IV. ABSOLUTE VS. SCALED POSITIONAL INFORMATION
The concentrations of morphogens provide cells with information about their position in the embryo.This "positional information" [30] can be measured in bits if we have access to data on the mean and variability of spatial profiles for the concentration of the relevant molecules [47,55].The local expression levels of individual gaps genes convey roughly two bits of information about position, twice what is possible in a model of on/off expression domains.Taken together all four gap genes provide ∼ 4.2 bits, sufficient to specify positions with ∼ 1% accuracy along the anterior-posterior axis, as seen above.However, these earlier analyses assumed, implicitly, that information is about the fractional or scaled position.Is this correct?
The key to separating information about scaled vs. absolute position is to compare the variance in morphogen concentrations at a scaled position x s depending on whether we constrain the length of the embryo (Appendix B).Qualitatively, if there is perfect scaling then knowing the length would not add any information with which to predict the morphogen concentration.Since information is mutual this would mean that all the available information is about the scaled position.To test this quantita- The extra information ∆I that Hb expression levels carry about absolute rather than scaled position, defined by Eq. ( 10) and evaluated from Eq. ( 14).Estimates are based on random choices of N embryos out of the full experimental ensemble (points; circles show means with standard deviations), and the extrapolation Nem → ∞ follows the methods of Appendix E (red line).The result is ∆I = 0.00 ± 0.008 bits (red circle with error bar).(B) The extra information ∆I conveyed by all four gap genes together, defined as in (A) by Eq. ( 10) but now evaluated using Eq. ( 15).Symbols as in (A); the result is ∆I = 0.038 ± 0.039 bits.Error bars are larger because we are analyzing a multidimensional code, but there still is no significant difference from ∆I = 0.
tively in the context of the gap genes, we have assembled data on N em = 301 embryos, in each of which we have reliable simultaneous measurements on the spatial profiles of expression in all four gap genes, as described in Appendix D.
Figure 3A shows the spatial profile of Hb as a function of scaled position along the anterior-posterior axis.At each scaled position x s = x/L we can visualize the distribution of expression levels, which is well approximated as a Gaussian (Appendix D and [47]).We can then ask if this distribution changes when we look only at embryos in a narrow range of lengths L, and the answer is no (qualitatively; Fig. 3B).Quantitatively we want to estimate the difference in entropy between these two distributions, averaged over x s and L, which will give us the deviation from scaling ∆I in Eq. ( 10), as explained in Appendix B. The calculation of the entropy simplifies in the Gaussian approximation, depending just on the variances as in Eq. (B19), where σ 2 g (x s |L) is the variance in concentration at scaled position x s across embryos of length L and σ 2 g (x s ) is the same variance computed across all embryos.
Applying Eq. ( 14) requires estimating the relevant variances and also making bins along the x s and L axes.For the scaled position we choose bins of size ∆x s = 0.01, consistent with the precision that we see in Figs. 1 and 2. To sample the range of embryo lengths we use N bins = 5, 10, 15, or 20 adaptive bins, and find the same results in all cases (Appendix E).As is well known, estimates of entropy or information are subject to systematic errors [38,56].In the present case, if we substitute estimates of the variances into Eq.( 14), we find a nonzero result for ∆I.But suppose we include different numbers of embryos in our analysis.In that case, we see that our estimate of ∆I depends on 1/N em as expected theoretically [38,56], and having seen this predicted dependence we can extrapolate N em → ∞.In particular, if we shuffle the data so that the true ∆I = 0, then our estimation procedure returns a random number with zero mean and standard deviation equal to our quoted error bar, demonstrating that we have control over the systematic errors.These now standard analysis methods are explained more fully in Appendix E.
Results of this analysis for Hb are shown in Fig. 4A.Using all N em = 301 embryos in our data set produces a very small estimate of ∆I, but even this is exaggerated by systematic errors as we see by changing N em .Our best estimate extrapolates to zero as N em → ∞, with an error bar smaller than 0.01 bits.When we repeat the same analyses for each of the other gap genes (i.e., Gt, Kni, and Kr), we get the same result (Appendix E).
We can generalize this analysis to consider all four gap genes simultaneously.Now the role of the variance in Eq. ( 14) is played by the covariance matrix of the fluctuations, as in Eq. (B23):  10) and evaluated from Eq. ( 14).Symbols as in Fig. 3, but the extrapolation now leads to a significantly nonzero value of ∆I = 0.1 ± 0.02 bits.Data from [49].
L, and Σ(x s |L) is the covariance computed across all embryos.Because we are looking at higher dimensional variations the impact of the finiteness of our data set is larger, but again we see the predicted dependence on 1/N em and can extrapolate to give ∆I = 0.038±0.039bits (Fig. 4B).
Once again this is consistent with ∆I = 0: there is no significant evidence for encoding of information about absolute, as opposed to scaled position.
Although the number of bits has meaning, it is useful to express the deviation from perfect scaling as a fraction of the information available about scaled position [47,55], Thus the patterns of gap gene expression scale with 1% accuracy, not just at discrete positional markers but across the entire range of graded spatial variations.

V. MATERNAL INPUTS DO NOT SCALE
Having observed scaling in the pair rule stripe positions and followed this back to the gap genes, it is natural to ask if we can go one step further and trace the scaling behavior of the patterning system to the maternal inputs.Of the three maternal inputs that drive patterning along the anterior-posterior axis of the fly embryo, much attention has been given to Bicoid (Bcd) [24,25].The protein is present at high concentrations in the anterior, and there is a nearly exponential decay of concentration with distance toward the posterior; one can monitor the dynamics of Bicoid protein concentrations quantitatively in live embryos using GFP-fusions [34].
Comparison across closely related species of flies shows that the length scale of this exponential decay varies in proportion to the mean length of the embryo [57].Insertion of bicoid genes from other species into Drosophila melanogaster produces protein concentration profiles with length scales appropriate to the host, but these are not sufficient to rescue the embryo from deletion of the native Bcd [58].These results emphasize the subtlety of comparison across species and the impact of genetic variations, leading us to re-examine the behavior of Bcd profiles across a large number of live embryos drawn from the same inbred laboratory strain used in the analysis of gap and pair rule genes.
Figure 5 analyzes Bcd profiles from N em = 582 live embryos [49].Measurements are taken during a small temporal window in nuclear cycle fourteen [34], and the only normalization (as with the gap genes) is to subtract a common background level from all the embryos and set the highest mean concentration to one.When we group the embryos into eight classes based on their length L, we see that the average concentration profiles in all groups are the same when plotted vs. absolute position, except for small effects at the posterior pole (Fig. 5A).If we plot vs. scaled position the different groups of embryos separate significantly (Fig. 5B), providing direct evidence against scaling.We make this precise using the same information theoretic approach as above and now find a significant nonzero value of ∆I = 0.1±0.02bits (Fig. 5C).This may seem like a small number, but this is related to the ∼ 4% scale of variations in embryo length.We conclude that the maternal inputs do not scale, in agreement with earlier suggestions [18].
We emphasize that the absence of scaling in the maternal inputs should not be interpreted as a form of noise.Indeed, absolute concentrations of Bcd protein are highly reproducible across embryos and this can be traced to highly reproducible numbers of mRNA molecules [49,59,60].Instead, we should think of the maternal inputs as a nearly deterministic response to the boundary conditions in the embryo, which also have a direct impact on the gap genes; see Eqs. (19,20) below.

VI. SCALING AND ZERO MODES
The results here strongly support the view that patterns of gap gene expression are genuinely scale invariant and that this is an emergent property of the gap gene network.Here we take the precise mathematical notion of scale invariance literally and explore its implications.While we do not pretend to have a detailed model, it is useful to have in mind a class of models for how patterns might form.As a caution we recall Turing's introductory remarks [14]: "This model will be a simplification and an idealization, and consequently a falsification." If we focus on variations just along the anteriorposterior axis x, and ignore the discreteness of nuclei, then the concentration g i of protein encoded by gene i plausibly obeys an equation of the form Here D i is the diffusion constant of species i, R i is the maximum rate at which these proteins can be synthesized, τ i is their lifetime before being degraded, and F i (g) describes all the potentially complex interactions by which all the proteins regulate the expression of gene i.We assume that the mRNA and protein dynamics have separate time scales so that one is effectively slaved to the other and we can write only one variable for each gene.Further, we neglect time scales that might arise in the process of regulation itself, such as switching between different regulatory states, so that F i (g) is an instantaneous function of the relevant concentrations.These assumptions are quite conventional, and other than this what we have written is very general.For example, the function F i (g) could describe both activating and repressive interactions, and these interactions could be combinatorial.These equations include as special cases Turing's original models [14] and their intellectual descendants [61,62].The maximum steady state concentration of each protein is R i τ i , and we can choose units in which this is equal to one, as with the normalized profiles of gap gene expression in Fig. 2A-D.For simplicity we will assume that all the decay times are the same, τ i = τ , although this is not essential for what follows; finally, we choose units of time such that τ = 1.Then we have where the length scale λ i = √ D i τ .This describes an autonomous network, which is not quite realistic for the gap genes-which are driven by maternal inputs-but should be sufficient to draw qualitative conclusions about the implications of scale invariance.
The length of the embryo appears not in the dynamical equations but in the boundary conditions.For most proteins, there can be no diffusive flux of molecules into or out of the ends of the embryo, so that The situation for maternal inputs is different; as an example, in making the egg the mother places mRNA for the protein Bicoid (Bcd) at the anterior end (x = 0), and this is translated continuously, so that where T Bcd is the rate of translation in appropriate units.Let us imagine that the final pattern we observe is in steady state, so that where the notation reminds us that length dependence can arise once we impose the boundary conditions.If we have true scale invariance as in Eq. ( 1) then if we make a small change in the length of the embryo, so that L → L + δL, the expression levels should change as but Eq. ( 21) still must be true.This requires that This seemingly abstract condition has a direct implication for the dynamics of the network.Suppose that the system is close to its steady state so that we can write and δg is small.Then we can linearize the dynamics in Eq. ( 18) to yield We recognize the term in brackets as the same one that appears in Eq. ( 24).To understand this connection it is useful to think of all possible spatial patterns of gene expression as points in a high dimensional space.Concretely we can write where the functions {ϕ µ i (x)} are the spatial "modes" of expression and the set {a µ } provides the coordinates of one expression profile in this multidimensional space.The number of modes is the number of genes times the number of independent points along the x axis, e.g. the number of rows of cells; for the gap genes the result is that the space has a dimensionality d > 300.We can choose these modes as eigenfunctions of the operator that appears in both Eqs. ( 24) and (26), where λ µ ≥ 0 means that the steady state is stable.Then so long as the deviations from the steady state are small, the dynamics of the network are simple in this coordinate system, Through Eq. ( 24) we see that perfect scale invariance implies a "zero mode," a particular mode of gene expression associated with eigenvalue λ µ = 0. Importantly this is not the steady state pattern itself, but an additional mode.
The existence of a zero mode has several implications: • Most literally, one component in the spatial pattern of gene expression will relax very slowly to its steady state, much more slowly than other components.Formally the relaxation should be as a power of time rather than exponential.
• The dynamics describe a "restoring force" that pulls the patterns of gene expression toward their steady state values; the eigenvalues are the spring constants associated with these restoring forces.Along the zero mode, there is no (linear) restoring force, and in the presence of any finite noise, the fluctuations along this mode will be very large compared with other modes.
• Along directions with nonzero λ µ the fluctuations in g will be approximately Gaussian so long as they remain small, as we see for the gap genes.But along the zero mode, there should be some deviation from Gaussian behavior.
There is evidence that the spatial patterns of gap gene expression can be compressed into a lower dimensional space, consistent with the idea that a zero mode dominates the dynamics [63].The (4 × 4) covariance matrix of fluctuations in gap gene expression is dominated by a single mode at almost all locations along the anteriorposterior axis, this large variance mode relaxes ∼ 10× more slowly than the lower variance modes, and one can even see hints of non-Gaussian behavior [26].
The existence of a zero mode is a statement about the linearized dynamics.If the absence of a linear restoring force continues for finite deviations from the steady state then there is a line of attracting spatial patterns rather than a single stable pattern.Different points along this line are the patterns appropriate to embryos of different lengths, and the final pattern is selected by boundary conditions.Line attractors have long been discussed for neural networks [64].It has been noted that models of the gap gene network might support such line attractors [28], and there are also suggestions that internal dynamics of the network can generate approximate scaling [29].The observation of nearly perfect scale invariance in the real network leads us to a much sharper version of these ideas.

VII. DISCUSSION
Scale invariance is an appealing concept.It quantifies the intuition that organisms are built from parts that are in proportion to one another, independent of an individual organism's overall size.There is a long history of searching for such scaling not just in adult organisms but at early stages of development, and the fruit fly Drosophila melanogaster has been a particular target for these studies [19-21, 29, 40].If we compare related species of flies we can see spatial patterns of gene expression that scale, on average, across 10× changes in embryo length [57,58], and similar results are obtained within a single species but with artificial selection for length variation [22].It has always been less clear whether scaling occurs without such large genetic variations, across the natural length variations in a single species.
We have explored scaling across many embryos from a quasi-inbred laboratory stock, minimizing genetic variation.Across this ensemble, we see length fluctuations with a standard deviation of ±4% but embryos in the tails of the distribution have lengths ±10% from the mean (Fig. A1).Following previous work, we measured the positions of discrete markers-such as the CF position, the peaks of pair-rule stripes, and the boundaries of gap gene domains-and found precise scaling of the absolute positions with embryo length.This is consistent with previous results, but what is new is the precision that we observe: markers are at positions that are scaled relative to the embryo length with an accuracy of ∼ 1% across the full extent of the anterior-posterior axis.This observed precision excludes a broad class of models that combine information from both ends of the embryo without explicit scaling [39][40][41][42].
There remains a gap between the positioning of discrete markers and the fuller notion of scale invariance.The gap gets smaller as we track more markers across a wider range of positions, but it would be attractive to address scale invariance directly.We have introduced an information theoretic approach that analyzes the full, graded spatial profiles of gene expression and measures similarity in the natural units provided by the intrinsic noise levels of these profiles.Concretely, we introduce a decomposition of the information that morphogen concentrations provide about position into a component about scaled position and a deviation from scaling.Applied to the gap genes in the early fly embryo, the result is clear: the deviation from scaling is less than one percent of the total positional information.It is perhaps surprising that we can make such a precise statement about the functional output of a complex network.
In contrast to the results for the gap genes and the pair-rule genes, at least one of the maternal inputs, Bicoid, does not exhibit scaling.We can see this "by eye," simply plotting profiles vs. absolute or scaled position, and these impressions are quantified by the same information theoretic approaches used to demonstrate scaling in the gap genes.Error bars again are in the range of ∼ 0.01 bits, but the deviation from scaling now is ∼ 10× as large.The conclusion is that near-perfect scale invariance is an emergent property of the gap gene network.
If we take scale invariance as a precise mathematical statement then the dynamics of the underlying genetic network must have a zero mode.This is equivalent to saying that the dynamics do not have a single attractor, but rather a line of attractors as in models for short-term memory in neural networks [64].Then position along this line is chosen by the boundary conditions and hence the length of the embryo.A zero mode would provide connections among several otherwise disparate observations on the gap genes.
Finally, recent experiments on mammalian pseudoembryos suggest that scale invariance may be a more universal feature of genetic networks underlying developmental pattern formation [65].In these self-organizing cell aggregates derived from stem cells, scale invariance emerges without fixed boundary conditions, but rather with boundaries that move as the aggregate grows.The existence of a zero mode in the regulatory network becomes even more attractive as a general mechanism for scaling.

I({g
Substituting, we can write where the three components are ) dx s dg P (g|x s ) log 2 [P (g|x s )] (B9) We can identify the three terms: First, I 1 is the information that the concentrations of morphogens provide about the scaled position, Second, I 2 is the entropy of the distribution of concentrations at a particular value of scaled position x s , averaged over this position, where S[Q] denotes the entropy of the distribution Q.
Finally, I 3 is the negative of the entropy of the same distribution but restricted to embryos of length L, and then averaged also over L, Comparing with Eq. (10) we see that the deviation from scaling can be written as the difference between two entropies, suitably averaged: This has a very simple interpretation: There is a deviation from scaling if specifying the length of the embryo reduces the entropy of fluctuations in morphogen concentration at a given scaled position.These general expressions simplify enormously in the case where we have only a single morphogen and the conditional distributions are Gaussian.In this case where ⟨g(x s )⟩ is the mean and σ 2 g (x s ) is the variance of g at scaled positions x s .Importantly the entropy of a Gaussian distribution does not depend on the mean, and we have [37,38] Thus we find the deviation from scaling is where σ 2 g (x s |L) is the variance in concentration at scaled position x s across embryos of length L. In other words, there is a deviation from scaling if the variance in morphogen concentration is reduced by knowing the length of the embryo.
This result can be generalized to multiple morphogens if we stay in the Gaussian approximation.Now with d genes, we have where Σ(x s ) is the covariance matrix of fluctuations in concentration at scaled position x s and ||Σ(x s )|| is the determinant of this matrix.Following the same logic as in the case of one gene we have (B23) Even if P L (g|x s ) is perfectly Gaussian, averaging over L could generate non-Gaussian behavior in P (g|x s ).We are neglecting this here, but since we find that ∆I is very small, and for the gap genes consistent with ∆I = 0, both the conditional and averaged distributions are very nearly Gaussian, as seen in Fig. 4B.

Appendix C: Cephalic furrow and scale invariance
Upon the onset of gastrulation (i.e., three hours after fertilization), the cephalic furrow (CF) emerges as the first macroscopic morphological feature along the anterior-posterior axis in the developing fly embryo.It results from collective cell movement that can be seen using bright-field microscopy.There are hints in early experiments that this marker is positioned very precisely [25].Modern experiments show that, as a fraction of the embryo length L, CF position x CF is reproducible to nearly 1% accuracy [49].
When we plot CF position in absolute units as a function of embryo length we observe a linear relationship with zero intercept (Fig. A2A).The slope of this fit is well within 1% of the mean scaled position ⟨f CF ⟩.More generally, all the discrete positional markers that we track (CF, pair-rule stripes, gap boundaries) have absolute positions that vary linearly with embryo length; the intercepts of the best-fit linear relations are zero; the slopes match the mean scaled positions of the markers (Fig. A2B) as predicted by scale invariance [Eq.( 3)]; and the precision of this match is better than 1%.

Appendix D: Aspects of the gene expression data
We analyze gene expression patterns for the pair-rule genes, the gap genes, and the maternal input Bicoid.In each case, the concentration of the protein is inferred from the intensity of a fluorescence signal.In each case images are collected by focusing on the midsagittal plane, the extent of the embryo is defined by thresholding the fluorescence intensity, and to avoid geometric distortions we avoid the 5% of the embryo at both the anterior and posterior poles.Fluorescence intensities are averaged over a small area, as shown in the inset to Fig. 1B, sliding along the dorsal rim of the embryo.In live embryos, we can keep track of time during nc14 directly, while in fixed embryos we use the progress of the cellularization as a clock with precision ∼ 1 min.
For each gene i we measure an intensity I α i (x s ) as a function of scaled position in embryo α.In each experiment, we normalize by assuming that the minimum mean concentration is zero and we choose units such that the maximum mean concentration is one.This defines where the background is and the scale factor is Importantly there is no freedom to normalize the profiles measured in individual embryos, which would distort our estimates of noise and variability [59].Data on three pair-rule genes-Eve, Prd, and Runare taken from Ref. [33].We fit the sum of seven Gaussians to each profile, identifying stripe positions with the centers of the Gaussians.We have also used the average peak profiles as templates [23] and made more restricted fits to small segments of the peaks [27]; results are the same with all three methods.Corrections for the drift of peak position vs. time in nc14 were made as described in the main text.
Simultaneous measurements on all four gap genes also were drawn from experiments described in Ref. [33].Because the analyses done here are so demanding of data, we tried to merge data from as many independent experimental sessions as possible.Most quantities are extremely reproducible from session to session, but in a handful of sessions, we found anomalously large variations in background fluorescence across the individual embryos.Concretely, if we measure the variance of expression levels for the individual genes, we typically find that σ 2 g (x s ) ≪ 10 −3 in regions of x s with near zero mean (Fig. A3).In a few sessions, these fluctuations in background are much larger, and these sessions are excluded;  more precisely, since all genes are measured simultaneously, excess background variance in one gene is sufficient to exclude those data.This leaves five independent sessions with a total of N em = 301 embryos which we pool into one data set for all further analyses.
For the analysis of gap gene expression boundaries, we mark the points that are half-maximal along the sharp slopes, as indicated in Fig. 2. For the weak peak of Kni expression near x s = 0.33 we fitted a Gaussian profile and took the marker as the center of the Gaussian.
Gap gene profiles vary slowly but significantly throughout nc14.If we don't treat this variation explicitly it can be confused for noise, resulting in a substantial overestimate of the variances and entropies.To separate the temporal dynamics of the gap genes from their noise level, we follow Ref. [32] and detrend the variations at each position, generalizing the treatment of the stripe positions in Fig. 1D.The alternative is to focus only on a small window of time [33], but this limits the size of the data set we can use.
Another systematic source of variation is the dependence of gap gene profiles on the dorsal-ventral coordinate [32].Previous work thus has been very strict in analyzing embryos with narrowly defined orientations.To expand our data set we are less strict, but this is problematic for the Kni profiles in the range 0.15 < x s < 0.45, which contains a small peak.When analyzing Kni alone, or all four gap genes together, we exclude this region.The alternative is to analyze the other three genes together across the full length of the anterior-posterior axis; results for ∆I are the same.An important assumption for our analysis is that the distribution of gene expression at a given anteriorposterior position is Gaussian, as shown previously [47].For completeness, we reproduce this result for our larger data set.In a single embryo α we observe a gene expression level g α (x s ) at scaled position x s .We compute the mean and standard deviation across all the embryos in our ensemble and define normalized deviations or z-scores We pool across all α = 1, 2, • • • , N em embryos and across all positions x s to estimate the probability density P (∆).
Results are in Fig. A4 for each of the four gap genes.Finally, measurements of Bicoid concentration are taken from fluorescent imaging of live embryos expressions a Bicoid-GFP fusion [49]; we consider only strain 2XA, which has the Bcd dosage closest to that of wildtype flies.With live measurements we can choose a time window, in this case t = 16 ± 2 min after the start of nc14, avoiding any temporal detrending while still including N em = 582 embryos.Some measurements with missing data points along the length of the embryo were excluded from this set.
Appendix E: Estimating ∆I from limited data Entropy and information depend on probability distributions, not just on their moments, and thus are especially difficult to estimate from a finite sample of data.Further, the entropy is a nonlinear function of the probabilities, and so random errors in probability become systematic errors in our estimates of information.This problem was appreciated in the very first attempts to use information theory in the analysis of experiments on biological systems [56].In the subsequent decades, many approaches have been developed, driven especially by the analysis of neural coding.The approach we take here follows the discussion in Appendix A.8 of Ref. [38].Rather than just saying that we follow established methods, we repeat some details that can be found in other contexts in hopes that our presentation thus will be more accessible.
If we estimate an information-theoretic quantity such as ∆I in Eq. (B14) based on data from measurements in N em independent embryos, then with any simple estimation procedure our estimates will be biased: Here ∆I ∞ is the true value of ∆I which we would observe if we could collect an infinite number of samples.The notation reminds us that if we make bins along some continuous axis, then the size of the corrections at finite N em depend on the number of bins N bins .With more bins the corrections are larger, which means that a naive estimate with a fixed number of embryos will depend on the bin size.The hope is that we can find a regime in which the extrapolated ∆I ∞ is independent of N bins .
It is important that Eq. (E1) is not just a guess, but a prediction that can be derived theoretically.Theory also gives values for the coefficients A and B, but these depend on details such as the independence of samples; the form is more general.This suggests a strategy in which we vary the number of embryos that we include in our analysis and look for the predicted systematic dependence on N em .If we can see this behavior then we can feel confident in fitting to Eq. (E1) and extracting an estimate ∆I ∞ [38].
This estimation procedure is illustrated by Fig. 4 in the main text and by Fig. A5.When we vary the number of embryos that we include in our analysis, we can choose at random from the total number available, so we have a path also to estimating error bars (below).In Fig. 4A we analyze ∆I for the spatial profiles of Hb expression using the Gaussian approximation of Eq. (B19).We have to make estimates of the variance as a function of the scaled coordinate x s and the embryo length L. As explained in the main text, we choose bins of ∆x s = 0.01, consistent with the observed precision of the pair-rule stripes and with earlier work [33,47].Along the L axis we use adaptive bins, that is bins with boundaries chosen so that the number of embryos in each bin is as nearly equal as possible; these bins are chosen based on the full experimental ensembles, and not readjusted as we choose smaller samples at random.
In Figure 4A we have chosen N bins = 5 adaptive bins along the L axis, and we see a clean depending of ∆I on N em as predicted in Eq. (E1).The dependence on N em is dominated by the linear term A, although some curvature is detectable.The quality of the fit to Eq. (E1) is very good, and the extrapolation is to ∆I ∞ = 0 with a precision of better than 0.01 bits.
In Figure A5 we see how the plot of ∆I vs. N em changes as we vary N bins = 5, 10, 15, 20.With more bins, there are fewer embryos in each bin, which means that the minimum number of embryos that we need for a meaningful analysis is larger.Increasing the number of bins might reveal otherwise hidden information, but also increases the size of systematic errors.Comparing across the panels in Fig. A5 we see that at fixed N em the apparent ∆I increases with N bins , and if we didn't explore the full dependence on N em we might be tempted to conclude that we are indeed revealing extra information.But this is not the case, since the plots at all values of N bins extrapolate to zero within error bars.
One useful test of these extrapolation methods is to be sure that we arrive at zero information in those cases where we know that the answer must be zero.As an example, if we randomly permute or shuffle the data we can break correlations that lead to nonzero information.In this case, if we randomly reassign lengths L to the embryos, then we must have ∆I = 0.This is illustrated in Fig. A6, using the example of Bcd.Here the real data extrapolate to a nonzero value of ∆I (Fig. 5), and when we shuffle we still see significantly nonzero values at N em ∼ 100.But using the whole N em dependence we can see this extrapolates smoothly to zero, as it should.
An essential part of this analysis is the estimation of error bars.For reasonably large N em the systematic and random errors are additive, and the variance of random errors scales ∝ 1/N em as usual.This means that if we compute the variance of ∆I across random halves of the data, and divide by two, we should have a good estimate of the variance in ∆I based on all of the data.If this error bar σ ∆I is correct, and our extrapolation is consistent, then when the true ∆I is zero, as with shuffled data, we can form a z-score, z = ∆ ∞ /σ ∆I , and z should be Gaussian with zero mean and unit variance.We can test this because the extrapolation procedure involves taking random subsamples of the data, each of which generates a slightly different value of ∆I ∞ .Figure A7 shows the distribution P (z) obtained from a shuffled version of the Hb data, illustrating a good match to the expected Gaussian distribution; the deviation is a bias toward smaller z, suggesting that our estimates of the error bars may be a bit conservative.Now that we have control over both the systematic and

FIG. 1 .
FIG. 1. Precise scaling of pair-rule stripes in the Drosophila embryo.(A) Bright-field image overlaid with fluorescent antibody staining for Eve protein (fuschia), focusing on the mid-sagittal plane with the dorsal side up; scalebar is 100 µm.(B) Expression of Eve in the second half of nuclear cycle fourteen (nc14).Solid line is the mean, and shaded region is the standard deviation across Nem = 108 embryos in a time window between 30 and 60 min from the start of nc14.Inset shows a single nucleus with a white square (width 0.01L) used to average intensities.(C) Eve expression profiles as a function of relative position along the body axis for 12 time bins during nc14, as indicated by color.(D) Linear dynamics of Eve peak positions during nc14, fit to Eq. (11).(E) Absolute positions of Eve peaks measured from the anterior pole referred to t0 = 45 min, as in Eq. (12), plotted vs. embryo length.(F) Standard deviation of scaled stripe positions as a function of mean position for three pair-rule genes, and for the cephalic furrow (CF, see Appendix C).Error bars are standard deviations from bootstrapping.Black curves with red shading (bootstrapped errors) are estimates of precision based on anchoring in Eqs.(5-7), and d is the spacing between neighboring cells.

FIG. 2 .
FIG. 2. Precise scaling of gap gene expression boundaries.Expression profiles of (A) Hunchback (Hb), (B) Giant (Gt), (C) Knirps (Kni), and (D) Krüppel (Kr), based on immunofluorescent staining (Appendix D).Means (solid lines) and standard deviations (shading) across embryos aligned by scaled position xs.Vertical lines indicate the mean positions of expression boundaries as well as a small peak in Kni.(E) Absolute position of all gap gene boundaries as a function of the embryo length.Dashed black line indicates the position of the posterior of the embryo.Boundary positions are time-corrected to t0 = 45 min, as with the stripe positions in Fig. 1D.(F) Standard deviation of scaled boundary positions as a function of mean position for all 13 markers.Error bars are standard deviations from bootstrapping.Black curves with red shading (bootstrapped errors) are estimates of precision based on anchoring in Eqs.(5-7), and d is the spacing between neighboring cells.Horizontal dashed lines denote the distance d and half-distance d/2, between neighboring nuclei.Dotted gray line indicates 1% precision.

FIG. 3 .FIG. 4 .
FIG. 3. Expression of Hb in scaled coordinates.(A)Mean concentration of Hb, ⟨g Hb (xs)⟩, vs scaled position (solid line, as in Fig.2A) and the conditional distribution P (g Hb |xs) around this mean (shading).Intensity bin size is 0.05 maximum ⟨g Hb ⟩. (B) A slice through the conditional distribution at xs = 0.47 (dashed black lines) compared with distributions estimated from embryos in narrow bins of length, PL(g Hb |xs).Lengths were binned in 5 bins with an equal number of embryos in each, such that each bin contains about 60 embryos with variations in L of less than 1%.Mean lengths in each bin are indicated at the upper right of each panel.Probability distributions of g Hb are estimated using a kernel density estimator with a Gaussian kernel that has width δg = 0.07 × maxx s ⟨g Hb (xs)⟩.

FIG. 5 .
FIG. 5.The maternal input Bicoid does not scale.(A) Measurements of Bcd concentration in Nem = 582 live embryos are grouped into eight classes by embryo length L and averaged.There is only one global normalization, so this shows that absolute concentrations have the same dependence on absolute position x across all classes.(B) The same data plotted vs. scaled position xs = x/L.Profiles separate, providing evidence against scaling.(C) Extra information ∆I that Bcd concentration provides about absolute vs. scaled position, defined by Eq. (10) and evaluated from Eq. (14).Symbols as in Fig.3, but the extrapolation now leads to a significantly nonzero value of ∆I = 0.1 ± 0.02 bits.Data from[49].
FIG. A2.Cephalic furrow and the proportionality of scaling.(A)The absolute position of the cephalic furrow measured in live emrbyos[49] is proportional to the embryo length.Red line is the best fit with 95% confidence intervals shown as dashed lines.The entire fit is shown in the inset, emphasizing that the intercept is zero.(B) Slopes of absolute position vs. embryo length for multiple positional markers, each plotted vs. its mean scaled position.(C) Replotting of data in (B) to show that slopes and scaled positions are equal within 1%, as predicted for perfect scaling [Eq.(3)].
FIG. A4.Fluctuations in gap gene expression are approximately Gaussian.Distributions of z-scored fluctuations, as in Eq. (D4), are estimated for each individual gap gene, pooled across positions and embryos; error bars are standard deviations.Black curves are Gaussians with zero mean and unit variance.
FIG. A5.Consistent estimates of ∆I with varying N bins .(A) Repeats the results of Fig 4A on ∆I vs. Nem for Hb, analyzed with N bins = 5 adaptive bins along the L axis.(B) As in (A) with N bins = 10; (C) with N bins = 15; and (D) with N bins = 20.Systematic errors are large at fixed Nem and increasing N bins , but the extrapolation Nem → ∞ is within error bars of ∆I = 0 in each case.

FIG. A6 .
FIG. A6.Recovering ∆I = 0 in shuffled data.We permute the lengths L of the embryos at random and repeat the analysis of the Bcd profiles.While the real data extrapolate to nonzero ∆I (Fig5C), here we recover ∆I = 0 as expected.

< l a t e x i t s h a 1 _
FIG. A7.Distribution of errors in ∆I.The distribution of errors of ∆I was estimated by repeating 5000 times the entire procedure leading to Fig.3A, on shuffled versions of the Hb data, where we know the true value of ∆I is 0 bits.We calculate the z-score based on our estimates of ∆I and error bars σ∆I .Red histogram shows the probability distribution P (z) in bins of size ∆z = 0.1.Shaded area is the uncertainty estimated through bootstrapping; black line is a Gaussian distribution with ⟨z⟩ = 0 and ⟨z 2 ⟩ = 1.